Analysis of spin precession in binary black hole systems including 
quadrupole-monopole interaction 



OO 
O 

o 

(N 
©/)■ 

< 



7 

OX), 



> 

o 

(N 
00 

cn 

o 

OO 
O 



X 



Etienne Racine 

Maryland Center for Fundamental Physics, Department of Physics, 
University of Maryland, College Park, MD 20742 
(Dated: August 11, 2008) 

We analyze in detail the spin precession equations in binary black hole systems, when the tidal 
torque on a Kerr black hole due to quadrupole-monopole coupling is taken into account. We show 
that completing the precession equations with this term reveals the existence of a conserved quantity 
at 2PN order when averaging over orbital motion. This quantity allows one to solve the (orbit- 
averaged) precession equations exactly in the case of equal masses and arbitrary spins, neglecting 
radiation reaction. For unequal masses, an exact solution does not exist in closed form, but we 
are still able to derive accurate approximate analytic solutions. We also show how to incorporate 
radiation reaction effects into our analytic solutions adiabatically, and compare the results to solu- 
tions obtained numerically. For various configurations of the binary, the relative difference in the 
accumulated orbital phase computed using our analytic solutions versus a full numerical solution 
vary from ~ 0.3% to ~ 1.8% over ~ 80 — 140 orbital cycles accumulated while sweeping over the 
orbital frequency range ~ 20 — 300 Hz. This typically corresponds to a discrepancy of order ~ 5 — 6 
radians. While this may not be accurate enough for implementation in LIGO template banks, we 
still believe that our new solutions are potentially quite useful for comparing numerical relativity 
simulations of spinning binary black hole systems with post-Newtonian theory. They can also be 
used to gain more understanding of precession effects, with potential application to the gravitational 
recoil problem, and to provide semi-analytical templates for spinning, precessing binaries. 



I. INTRODUCTION AND SUMMARY 

Coalescing black hole binary systems are a class of 
sources of gravitational radiation that can potentially 
be detected by observatories such as LIGO [l(, VIRGO 
0, GEO and TAMA 4]. Data analysis algorithms 
searching for such sources are based on the method of 
matched filtering, which essentially consists of correlat- 
ing the detector output with a bank of (theoretical) tem- 
plate signals. For detection to be at all possible, these 
templates must be accurate to ~ 1 radian over hundreds 
to thousands of cycles. In the regime where the black 
holes are separated by a large enough distance, the post- 
Newtonian (PN) expansion is expected to be accurate 
enough to be useful in generating gravitational wave- 
form templates, provided one carries the expansion to 
high enough (3.5 PN) order (e.g. see [j| and references 
therein). 

While the distribution of spins for black holes in bina- 
ries is still a matter of ongoing research, it is generally 
thought very plausible that these spins may be large. It 
is therefore important to include in template banks wave- 
forms that account for black hole spins, as they leave sig- 
nificant imprints on the signal, e.g. strong modulation 
of the waveform amplitude and phase 0, U 11 • Ignoring 
effects of black hole spins could very well lead to missing 
such binaries when analyzing interferometer data. 

In order to predict the contribution of the spins to 
the gravitational waveform to the accuracy required by 
the data analysis process, the time evolution of the spins 
must be computed at least up to 2.5 post-Newtonian or- 
der (e.g. see 0, [l(| and references therein). However 
in this paper we restrict attention, for simplicity, to the 



precession equations accurate at 2PN order, i.e. includ- 
ing leading order spin-orbit and spin-spin terms. To that 
order the spin evolution equations, which will also be 
loosely referred to as "precession equations" throughout, 
take the functional form 



dt 



fo{t) J x Si )2 + /i(t)S 2 ,i x Si, a 

+ terms of higher PN order, (1-1) 



where Sq,2 are the spin 3-vectors, J is the system's total 
angular momentum, and fo,i(t) are known, implicit func- 
tions of coordinate time. In the literature the functions 
/o,i are generally computed assuming a model of spin- 
ning point-particles, i.e. the precession equations are de- 
rived from solving the PN expanded Einstein equations 
assuming a distributional stress-energy tensor that con- 
tains mass monopole and spin dipole degrees of freedom. 
However this procedure ignores the precession due to the 
coupling of the mass quadrupolc moment of each spin- 
ning hole to the tidal field of the companion 1 . Since the 
mass quadrupole moment Q of a Kerr black hole scales 
in order of magnitude as Q ~ S 2 /M, where M is here 
the black hole mass and S the magnitude of its angular 
momentum, the precession of the hole's spin due its mass 
quadrupole contributes at the same PN order as the spin- 
spin term in Eq. fll.ip . but is nevertheless often omitted 
in the literature, with the following notable exceptions. 



1 This is essentially the general relativistic version of the well- 
known phenomenon of the precession of the equinoxes on Earth. 
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Let us first point out the works of Barker 
and O'Connell [ll[ and Damour [l^], who consider 
quadrupole-monopole coupling in the precession prob- 
lem. Barker and O'Connell derive generic precession 
equations from a Hamiltonian point of view, while 
Damour focuses to the problem of spinning binary black 
holes, which allows a very compact Hamiltonian formu- 
lation at 2PN order. Next we note the contributions of 
Poisson [l3| and Pan et al. Q, who have looked at some 
contributions of this term in the context of gravitational 
wave data analysis. In Ref . [13[ Poisson argues that since 
the quadrupole-monopole precession contributes at the 
same PN order as spin-spin precession, the main features 
of the work of Apostolatos et al., who study extensively 
spin precession in binary black holes at leading order in 
post-Newtonian theory, should remain qualitatively un- 
changed by the addition of the quadrupole-monopole pre- 
cession term, as it is a sub-leading term. In Ref.[H, Pan 
et al. devote a subsection to a comparison of gravita- 
tional wave templates with and without the quadrupole- 
monopole precession, their conclusion being that the con- 
tribution of the quadrupole-monopole precession can be 
neglected. Lastly it is worth mentioning the works of 
Gergely, Keresztes and Mikoczi [l4], [H| , who consider the 
quadrupole-monopole contributions to both conservative 
and dissipative pieces of orbital dynamics of compact bi- 
naries, and in particular provide a complete solution to 
the binary's radial motion. 



the discovery of a new conserved quantity 2 . We then 
show in section IIVI that this conserved quantity allows 
one to solve the orbit-averaged precession equations ex- 
actly in the case of equal masses and arbitrary spins, 
provided one neglects radiation reaction. In the case of 
unequal masses, we derive an accurate perturbative an- 
alytic solution. These solutions allow one to identify all 
the relevant frequencies characterizing the precession and 
nutation motions of each spin relative to the total angular 
momentum axis. Finally in section [Vl we derive approx- 
imate analytic solutions to the precession equations tak- 
ing into account radiation damping, and compare these 
solutions to numerically integrated solutions. 

II. QUADRUPOLE-MONOPOLE 
INTERACTION 

In Newtonian gravity, it is well-known that a body pos- 
sessing a mass quadrupole moment (or any higher order 
mass multipole moment) placed into a prescribed grav- 
itational tidal field (e.g created by a companion object) 
experiences a torque, which leads to a non-trivial time 
evolution of its spin. This effect is, for example, what 
causes the precession of the Earth's spin axis (precession 
of the equinoxes). A simple computation of the torque 
Tj in Newtonian gravity, using a coordinate system (t, x) 
that is mass-centered on the body, yields 



However in light of the analyses of Poisson [l3( and Pan 
et al. , it is understandable that a complete study of the 
precession equations including the quadrupole-monopole 
contribution for black hole binaries has not yet been 
completed. The main purpose of this paper is to per- 
form this analysis carefully, being motivated by the fol- 
lowing two considerations. The first element is consis- 
tency: in black hole binaries the quadrupole-monopole 
term contributes at the same PN order than spin-spin 
coupling. Thus if spin-spin coupling is taken into ac- 
count in constructing template waveforms, then so should 
quadrupole-monopole coupling in principle. The second 
motivating element is essentially academic interest in an 
unsolved problem, coupled with the fact that sometimes 
order of magnitude estimates miss important facts and 
more careful investigations reveal hidden features in a 
problem. In the analysis we present below, such a fea- 
ture will be unravelled, namely the existence of a new 
conserved quantity in the orbit-averaged precession equa- 
tions, when the quadrupole-monopole precession is taken 
into account. 



n = J d 3 x p{x,t)e i: jkX3[-dk$ cx t{x 7 t)} 



d 3 x p(x, t)eijkXj 



oo 

d k Y i ^G <ili2 ... il> {t)x< h x i \..x i ^>, (2.1) 



1=0 



where the brackets < ... > denote the symmetric trace- 
free projection, p{x, t) is the mass density of the body, 
$ ext (£c,t) is the Newtonian potential created by an ex- 
ternal source, and where 3 



G 



<«ll2— tl 



> (t) = -di 1 di 2 ...d il § e xt{x,i) 



(2.2) 



The right-hand side of (|2.2p is clearly symmetric, and 
Laplace's equation ensures that is also trace-free on all 
indices. Using the Newtonian definition of mass multi- 
pole moments 



Our analysis is organized as follows. We first review 
in section [TT] the derivation of the spin precession term 
due to the mass quadrupole of a Kerr black hole. In 
section IIIII we show that including the mass quadrupole 
term in the orbit-averaged precession equations leads to 



2 This quantity is strictly conserved by the orbit- averaged preces- 
sion equations if radiation reaction effects are neglected. If radi- 
ation damping is taken into account, it evolves over the radiation 
reaction timescale. 

3 Note that the 1 = term does not contribute in 112.11 1 due 
to the spatial derivative dk, and thus we do not need defining 
G <ili2 ...i l >{t) for I = in (2JJ. 
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/dS 1 
<Pxp{x,t)x <il x i '...x il> , (2.3) -j± = ^ 



we obtain 



1 f 

e ijk^2 n G kii-h / d 3 X p(x, t)x j X <%1 i 
1=0 ' J 

e ijk}]^Gkix...ii / d 3 xp(x,t)Q 

1=0 ' ^ 
e ijk^ -n M jii-..i l Gki 1 ...i l , (2.4) 



i=0 



where the second line follows from the fact that 
x J x <n x 12 ...x ll> and x <3 x 11 x %2 ...x tl> differ by trace 
terms, which do not contribute since £ijkGkix...j...ii 
0. The form (|2.4p of the torque follows the syntax 
of Damour, Soffel and Xu 16]. When specializing to 
a binary system and restricting attention solely to the 
quadrupole-monopole interaction, we find 



dSj 



<om> M 



n <k n m> 



(2.5) 



where Qf jm> is the mass quadrupole of body 1, M 2 is 
the mass of body 2, n l is a unit vector pointing from 
body 1 to body 2 and r is orbital separation. While 
Eq. (|2.4[) has been revisited here in the context of New- 
tonian gravity, its validity in the case of a binary black 
hole system in the regime where the post-Newtonian ex- 
pansion is valid has been proven rigorously by means of 
a surface- integral method by Racine in [17]. For that 
situation, the mass multipole moments of each object 
are defined through their imprints in the far-field metric. 
Thus the mass quadrupole moment entering Eq. (|2.5p in 
the case of a spinning black hole is well-known and is 
given by (e.g. see Poisson [l3|) 



(2.6) 



where S\ — xiMfS[, Xi being the dimensionless magni- 
tude of the black hole spin, and S\ the unit vector along 
the spin of body 1. Substituting (|2.6[) into (12.51) we obtain 



dS 2 _ 1 

dt ~ r 3 



2 + A 
2q 



' 3q 

2+ y 



L N - S 2 + j^[fi ■ S Q ]n\ x Si, 

(2.8a) 

L N - Si + |f [h ■ So]n\ x S 2 , 
1VI2 I 



(2.8b) 

where q = Mi/M 2 , p = MiM 2 /(Mi + M 2 ), and where 



So 



Mi 
W 2 



(2.9) 



The vector So was originally introduced by Damour in 
the 2PN Hamiltonian for spinning binary black holes [13] • 
To close the system (neglecting radiation-reaction for the 
time being), one uses conservation of total angular mo- 
mentum 4 J = Ln + Si + S2 to derive the precession 
equation for the (Newtonian) orbital angular Ln, which 
is 



dL 



N 



1 



dt 









■I*) 





■iN 



3/z 

M 



(n • S )h x S 



(2.10) 

where S = Si + S 2 - One may also recognize the often- 
encountered combination 



1 3 
2S e ff = -S + -So- 



(2.11) 



III. A NEW CONSERVED QUANTITY 



From Eqs. (|2.8p . one can estimate the timescale for pre- 
cession of the spins r p to be r p ~ r 3 /\L^\ ~ (r/p)T OI t,, 
r or b being the orbital period. Thus it is customary to av- 
erage the precession equations (I2.8[) - ([2.10[) over an orbital 
period if one is interested in the secular evolution of the 
spins, i.e. neglecting small fluctuations occurring over 
the orbital timescale. This average over orbital motion 
is performed using (see e.g. Schnittman [191 ]) 



1 



1 



V a 3 (l-e 2 ) 3 / 2 ' 
,n l n J , 1 



2a 3 (l - e 2 ) 3 / 2 



(3.1a) 
(3.1b) 



f = £f(*-*)Ax* (2-7) 

for the quadrupole-monopole precession term. When 
adding this term to the precession equations typically 
found in the literature (see e.g. 0, El), the Ml preces- 
sion equations are 



4 Strictly speaking, the true conserved angular momentum at 2PN 
is expanded as J = L + S = Lpi + Lpn + I>so + I*2PN + S (sec 
e.g. [31 ). However we can neglect all post-Newtonian corrections 
to L, i.e. we take L = -Ljv, since these corrections introduce 
terms in the precession equations which are of higher PN order 
than the terms considered in this paper. 



4 



where a is the semi-major axis and e is the eccentric- 
ity. Defining d = aV 1 — e 2 , we then obtain the following 
system of evolution equations 



dSi 
~df 

dS 2 
dt 



dL 



N 



dt 

where 



1 

2d? 

1 

2d? 



2d 3 



Mi Mi 
' 3Mi 3^ 



L N + S 2 } x Si, 

(3.2a) 



M 2 M 2 



{s + 3 



1 + A A 
M 



A = 



1/JV ' Sq 



A + Si j x S 2 , 

(3.2b) 

So} x L N , (3.2c) 



(3.3) 



A key property of system (|3.2p is that it possesses a previ- 
ously unknown conserved quantity, whose identification 
was only possible by completing the precession equations 
with the quadrupole-monopole interaction. This con- 
served quantity turns to be A, defined in Eq. (|3.3p . This 
can be shown by a direct computation. First, clearly 



'N\ 



is conserved from (|3.2cl) . giving 



dX ( dS , „ dL 

n a \ LN -^r + s °- 



N 



dt 



cx \M 2 L N • (S 2 x Si) + MiL N • (Si x S 2 ) 

+/iS ■ (S x L N )} 
= [M 2 L N • (S 2 x Si) + M X L N • (Si x S 2 ) 

+(M 2 Si + M X S 2 ) • [(Si + S 2 ) x Ljv)]} 
= [m 2 [L n ■ (5 a x Sj) + Si ■ (S 2 x Ljv)] 

+M 1 [i JV • (Si x St) + S 2 • (Si x L N )]j : 



= 0. 
(3.4) 



It is important to note here that the conservation of A 
is a property of the orbit-averaged precession equations 
only. Repeating the same computation using (|2.8p and 
(|2TTD|) instead of (MD yields dX/dt ^ 0. This implies, for 
example, that the conservation of A (when averaging over 
orbital motion) cannot be deduced directly from the 2PN 
Hamiltonian of spinning binary black holes [l2j], unless 
the orbital degrees of freedom are averaged out. 

In addition, proof (|3.4p is only valid over the precession 
timescale. Over the radiation reaction timescale, the or- 
bital angular momentum evolution equation must be sup- 
plemented by the following dissipative term [restricting 
henceforth attention to orbits with negligible eccentricy, 
i.e. d-> a + 0(e 2 )} 



• _ 32 fiM 2 



(3.5) 



The corresponding evolution equation for A is then easily 
shown to be 



dX 32 fj,M 
~dt ~ ~5 dF 



r2 



■A = 7 A, 



(3.6) 



whose solution can formally be writt en as t he following 
integral over orbital frequency uj = V Ald~ 3 



X = Xq exp 



7(cj) 



djj 



(3.7) 



where Aq and loq are initial values. 



IV. ANALYTIC SOLUTIONS IN ABSENCE OF 
RADIATION REACTION 



The first step in solving the precession evolution equa- 
tions is making use of the definition of total angular mo- 
mentum J = Ln + Si + S 2 to eliminate Ln from (|3.2a|) 
and (|3.2bp . We then have 



dSi 
~df 



dS 2 
dt 



1 

2d? 



1 

2d? 



3M 2 
Mi 
3M 3/^ 
M" ~ M~i 
3Mi 

3[i 



~ Ah 
S 2 | x Si, 
3// 



-A 



J 



3M 



M 2 M 2 



M 2 
Si) x S 2 , 



(4.1a) 



(4.1b) 



where J is a constant vector. The usefulness of the con- 
served quantity A is now quite clear: it shows that the 
coefficients in the coupled evolution system for Si and 
S 2 are all constants. 

In the remainder of this section we solve (|4.ip for both 
equal mass and unequal mass cases and arbitrary spins, 
generalizing the analytical analysis of Apostolatos et al. 
[6j, valid for equal masses, or when one of the spins is 
dynamically negligible. 



A. Equal mass case 

When Mi = M 2 , the structure of system (I3.2[) simpli- 
fies greatly, allowing in fact an exact solution. Defining 
the phase variable i/j as 5 



We keep the integral over t' explicit here, in anticipation of our 
treatment of radiation reaction in the next section. 



5 



and the parameter a as 



Now S 2 , J • S and S • A are all constants of motion, 
which are thus known in terms of input initial conditions. 
(4 2) However J ■ A is not yet known. Hence before (|4.8j) can 
be solved explicitly, we must first solve for J ■ A. This is 
done directly by contracting (|4.8[) with J, giving 



4- A 



14- 3A 



the evolution of the spins is governed by 



(4.3) 



d 2 (J-A) , 2a2 



dtp 2 

The solution is immediate 



S (J ■ A) — a (S ■ A)(S • J). (4.9) 



dS 

— — = J x Si - aS 2 x Si, 



rfS 2 



j x S 2 — aSi x S 2 



(4.4a) 
(4.4b) 



By adding these two equations together, we see immedi- 
ately that the total spin S = Si + S 2 simply precesses 
about J, i.e. 



dS_ 

dtp 



= JxS. 



(4.5) 



The solution to differential equation (|4.5p is 



S = S l < + S£ cos tp + ( J x Sq) sin tp, (4.6) 

where Sq and Sq are the projections of the initial total 
spin So [not to be confused with the combination of spins 
defined in Eq. (|2.9p ; henceforth So will always refer to ini- 
tial total spin] in directions parallel and perpendicular to 
J respectively. Equation (|4.6|) is the first part of the ex- 
act solution for the equal mass case. We again emphasize 
that this solution corresponds to simple precession of the 
total spin about the total angular momentum axis, with 
precession frequency Q p = dtp/dt. 

Consider next the vector A = Si — S 2 . Its evolution 
is governed by 



dA 

~dip 



= J x A - 2aS 2 x Si 

= J x A - aS x A 
= n x A. 



(4.7) 



Taking one more derivative with respect to tp yields after 
some algebra 



d 2 A 

dtp 2 



(J - aS) 2 A = (J • A) (J - 2aS) + a 2 {S ■ A)S. 

(4.8) 



(J ■ A) = Acoaa + B sma + {S ■ A)(S ■ J), (4.10) 



where the phase variable a is defined below in Eq. (|4 . 19[) . 
To fix the integration constants A and B, we use the 
initial conditions and Eq. (|4.7p to obtain 



(J -A) = vjs(So ■ A ) + [(J x So) • ( A X S )] cos a 
+J- (A x S ) sine-, (4.11) 



where ^js = So ■ J ■ Equation (|4.11[) shows that the 
projection of the difference of the spins along the to- 
tal angular momentum axis oscillates at a frequency 
slower than the total spin precession frequency, since 
da/dt = (aS)Q p [cf. Eq. (|4~19)) below]. This implies that 
each individual spin Si ; 2 exhibits a nutation motion at 
frequency Q n = (aS)Q p . [Here we define nutation motion 
as non-trivial time evolution of the projection of a given 
spin vector along J.] Since the nutation frequency is 
proportional to a, it is entirely due to spin-spin coupling. 
Barker, Byrd and O'Connell [2fj] have computed spin nu- 
tation frequencies in generic binary systems. However it 
is difficult to compare their results to ours since they rely 
on a succession of approximations, along with a some- 
what implicit, time-dependent definition of the nutation 
frequency. In addition it is not entirely clear how their 
definition of nutation motion maps onto ours. Since our 
result is exact, up to the same PN order as considered by 
(20| . and provides a clean, constant value for the nuta- 
tion frequency along with a precise definition of nutation 
motion, we believe our result is an improvement over the 
work of [2(J. However our result is applicable only to 
binary systems of Kerr black holes. 

Substituting (|4.11j) into l|4.8p . one can then solve P~g]) 
using the retarded Green's function for the simple har- 
monic oscillator. However it turns out to be more con- 
venient to simply project (14. 8p along the remaining two 
independent directions, which we choose to be along Sq 
and J x Sq , in order to match the form of solution (|4.6p 
for the total spin. The remaining evolution equations are 
then 
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dip 



a (^-A) 



= -U 2 (S£- • A) + a 2 {S^) 2 {S Q • A ) cosiP - 2a(S L ) 2 (J • A) cos ip 

= -n 2 (S^ ■ A) + a(S^f(S ■ A )(aS - 2fx JS ) cost/- 

-a(^) a {[(Jx S )-(A x So)] [cos [(1 + aS)^] + cos [(1 - aS)^]] } 

+a(S L ) 2 { [J ' (ft x A )] [ sin [(1 + aS)V>] - sin [(1 - aS)iP] ] }, 



(4.12) 



^[(ix^).A] 



-n 2 [(J x ^o 1 

-|2 



The solutions to these forced harmonic oscillators are 



l (S£) 2 {S a ■ A ) sinV- - 2a(S L ) 2 (J • A) sin^ 

n 2 [(J x S£) ■ A] + a(S^) 2 (S ■ A )(aS - 2fi JS ) sin V 
a(S^-) 2 { [(J x S ) • (A x S )] [sin [(1 + aS)^) + sin [(1 - aS)if>]]} 

a(S^) 2 {[J ■ (So x A )] [cos [(1 + aS)i>] - cos [(1 - aS)ip]]\. 



(4.13) 



Sq ■ A = [S^ ■ A ] cosro + [S^ ■ (f2 x A )] sintu + 5(1 - /i 2 s )< (S ■ Ao)(cos-0 - cosro) 



--[(JxSo)-(AoxSo)] 



cos w — cos(tp — a) cos w — cos(^ + er) 



;[j-(S X A )] 



(1-Mjs) + 
(5+ sin zd — sin( - + a) 5— sin w — sin("0 — a) 



(1 + Ws) 



(1 - Mjs) 



(4.14) 



(J x Sq) ■ A = [(J x S^) ■ A ] costu + [(J x Sq 1 ) • (O x A )] sintn + S(l - /4sK (S • A )(sin ?/> - ft" 1 sintn) 



+ 2^ x So)-(A x So)] 



5- sin tu — siri(^ — a) 5 + sin w — sm(ip + a) 



-[J -(S x A )] 



(l-/ijs) (1 + Mjs) 

cos(ip + a) — cos nj cos("0 — <r) — cos w 



(1 + MJs) 



(1 ~ Mjs) 



(4.15) 



r 



where 

il = = (1 -2aSfi JS + a 2 S 2 ) 1/2 , (4.16) 

8± = n- 1 (l±aS), (4.17) 

and where the phase variables w and a are defined as 



— 



ndtp', 



aS dip' ■ 



(4.18) 



(4.19) 



Here we note that the components of the difference of 
the spins perpendicular to the total angular momentum 



axis contain terms which oscillate at four different fre- 
quencies. One of these frequencies is again f2 p , and the 
other three are close, but not equal to the total spin pre- 
cession frequency fl p , namely (1 — 2aSfijs + a 2 S 2 ) 1 ^ 2 Q p 
and (1 ± aS)fl p . Spin-spin coupling is again responsible 
for the appearance of these three additional precession 
frequencies. 

The complete vector A is then given by combining 

mum , (gnu) and dum to § ive 



(j-a)j+ J g ° x ; A L s, 



[(J x So 1 ) • A] 



JxS±. 



ll-^ 2 ° (1-^s)^ 2 

(4.20) 

The individual spins are recovered from Si. 2 = j(S± A). 
We emphasize once more an important property of this 
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solution: each individual spin exhibits nutation at fre- 
quency f2 n = (aS)Qp, which is typically much slower 
than the precession frequency of the total spin. Since 
the nutation frequency is directly proportional to a, this 
clearly shows that this effect is entirely induced by spin- 
spin coupling. The evolution of the components of the 
spins transverse to J is rather complicated, since it con- 
tains components oscillating at four different frequencies, 
which are ftp, (1 - 2olShjs + a 2 S 2 f/ 2 il p and fl p ± Q n . 
Again if one neglects spin-spin coupling, all four frequen- 
cies of the precession motion reduce to a single precession 
frequency, namely f2 p . 

Note also that in the case where, say, |S 2 | *C 
i.e. the spin of body 2 is dynamically negligible, then 
Si becomes the total spin S and clearly the vector A 
also becomes equal to the total spin. Obviously solution 
(|4.6|) for the total spin is still valid in that situation. It 
is straightforward to check that by setting Ao = So in 
(|4.20[) . one recovers A = S for all time. 



B. General case 

Unlike the equal mass case, the general case of arbi- 
trary spins and unequal masses unfortunately does not 
allow an exact solution in terms of elementary functions. 
We first prove this statement. The spin evolution equa- 
tions in the general case are 



where 



dSi 

Ihp 

dS 2 

d%p 



«1,2 = 



0i J x Si — aiS 2 x Si, 
(3 2 J x S 2 - a 2 Si x S 2 , 



Pi 



1 

\J 

4 + 3 



6 fM-fiX 



M h2 V 14 - 3A 



Ail,: 



Mi, 2 



7-§A 



(4.21a) 
(4.21b) 

(4.22) 
(4.23) 



Now consider the evolution equations satisfied by J ■ Si, 2 
and Si • S 2 . These are 



d(J-S^) 

dip 
djSi ■ S 2 ) 

dip 



±a h2 J • (Si x S 2 ), (4.24) 
/?_ J ■ (Si x S 2 ), (4.25) 



where /?_ = 0i — f) 2 . We see that all right-hand sides in 
the above equations are proportional to the same func- 
tion, namely J ■ (Si x S 2 ). Defining p as 



we then have 



J-Si, 2 = ±ai, 2 p+(J-S 1>2 )o, (4.27) 
Si-S 2 = (/3_)p+(Si-S 2 ) , (4.28) 

where the subscript denotes initial value. By taking 
a ip derivative of, say, Eq. (|4.25|) and using expressions 
(|4.27p - (|4.28p . one can show that p obeys the following 
differential equation 



d 2 p 

di\> 2 



-3P-aia 2 p 



ffL + 2/3_ J • C - + (C + ) 2 



-/3_ (J x Si) • (J x S 2 ) 



J x (Si x S 2 ) 



(4.29) 



where Ct ~ { a 2S\ ± aiS 2 )o- Thus p evolves as an an- 
harmonic oscillator with a cubic term in its potential. In 
the equal mass case (3- = and p reduces to a simple 
harmonic oscillator, thus allowing a solution in terms of 
trigonometric functions. In the general case one can solve 
(|4.29p in terms of an elliptic integral 6 , more specifically in 
terms of the Weierstrass elliptic function (see e.g. [2l|). 
However here, for sake of transparency and simplicity, we 
resort to a perturbative method to solve (|4.29p . based on 
the existence of a small parameter that will be identified 
shortly. Equation (|4.29[) is of the form 



d 2 p 

dip 2 



-C 2 p -Cip + Co, 



(4.30) 



where Cb,i j2 are all constants, which can be read off di- 
rectly from (|4.29p . Defining Si j2 = \S\ 2\ and the dimen- 
sionless variable u = p/SiS 2 , we have 



a u 
'chp 2 



-SiS 2 C 2 u 2 -Ciu4 
-B 2 u 2 - B lU + B Q . 



Co 
SiS 2 



(4.31) 



Now the key element to realize is that B 2 is a small num- 
ber, as it scales as 



Bo 



(aiSi)(a 2 S 2 ) 



M 2 Si S 2 
J 2 Mi M 2 



MiM 2 —j-xiX2 



M 2 
J 2 

V-J2X1X2, 



(4.32) 



p= [ J (Si x S 2 )dip', (4.26) 
Jo 



6 The author thanks Saul Teukolsky for pointing this out. 
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where (|4.22[) has been used and where X1.2 — S\,2 /M^ 2 
are dimensionless spins and rj — \ijM is the symmetric 
mass ratio. 

Let us first consider the case where 77 — > 0, i.e. the 
extreme mass ratio limit. We adopt the convention here 
that M 2 — > 0, so Mi — ► M. In that case the total angular 
momentum is dominated by the total spin since £jv ~ f}- 
Therefore we have J 2 — > M 4 x 2 which gives 



B 2 



VX2/X1 



0, 



(4.33) 



since xi is finite. 

When masses are comparable in the regime where post- 
Newtonian theory is applicable, the total angular mo- 
mentum is dominated by the orbital angular momentum, 
which yields J 2 ~ r] 2 M 3 r. This then yields 



□ X1X2 M X1X2 2 
B 2 ~ ~ v . 



v 



(4.34) 



Since X1X2/V is typically a number of order unity when 
masses are comparable, B 2 scales as a 1PN correction 
term compared to the numbers B\ and i?o, which are 
generically of order unity. 

Provided that u remains small enough, the anharmonic 
term in the equation of motion for u generically remains 
a small correction to the harmonic driving term for all 
time. We shall assume this is the case in what follows, 
and verify a posteriori that this assumption indeed works 
quite well. We next define the variable 



X = U — Aq , 



(4.35) 



where 



An = 



-B x + y/ B'{ + 4^52 
2B 2 



(4.36) 



is a constant chosen so that the evolution equation for x 
takes the form 



2 2 
-ex — w x, 



(4.37) 



where e — B 2 and 



W 



B 1 + 2B 2 A 



B\ 



4B B 2 



(4.38) 



We now solve (|4.37p by treating the anharmonic term as 
a small perturbation, with e being our expansion param- 
eter. Now it is well-known (see e.g. Goldstein [22]) that 
the corrections to the frequency of an oscillator with a 
cubic perturbation of order e in its potential scale as e 2 . 
Thus it is necessary to solve (|4.37p to order e 2 accuracy 
in order to capture this important effect. To that order 
the solution is 



+ 



x + v 1 X Q 

2x v , 2 



xo 



2v 2 



3w 2 



,xo(29a;o — 55w 2 ) 
144w 4 



cos($) 



V - £■ 



5^0(^0 



llv 2 ^ 



3w 2 



+e 



X V Vq( 

3ttr 



9w 4 



U4w 4 

sin(2$) 



sin(<i>) 



H) 



,xo(xl 



48w 



xpjxl + Av 2 ) 

3v 2 ) 9 uo(3a;n 
p-^ cos(3$) + e 2 v ; 



6w 2 



cos(2$) 



48w 4 



•sin(3$), (4.39) 



where vo = (dx/d$)o, with 



terms of initial conditions of the spins, the remaining 
quantities appearing in (|4. 39[) are 



Mw 2 x 2 + (x' ) 2 ] 
6w 6 



(4.40) 



where x' = (dx/dip)o. Since the components of each spin 
along J are proportional to x(ip), one can identify the 
fundamental frequency associated with the nutation mo- 
tion of each spin as (d&/dip)ttp. However in the general 
case, the nutation motion contains higher harmonics of 
its fundamental, as opposed to the equal mass case where 
the nutation motion contains exactly one frequency. In 



x = 



2B 2 



,5[w 2 x 2 + (x' ) 2 



6w 6 



-1/2 , 

^, (4.42) 
w 



x' = J (Si x 5 2 )o, 



(4.43) 
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w 



[Bl + 4B B 2 



1/4 



(4.44) A 1 - in (|4.49[) by their average 7 . The resulting system of 
equations for S and A 1 - is then of the form 



Bo = Co ■ [J x (A x S; 



2JJ0 



-04(JxS 1 )-(JxS 2 )}o, (4.45) 



Bt = 2 _ + 2/3_ J ■ Co + (Co + ) 2 , (4-46) 



B 2 = e = 3(0_)(a 1 S 1 )(a 2 S 2 ). (4.47) 

Now recall that the function x(ip) [or equivalently p{ip)} 
allows one to determine only three angles in the problem, 
namely the projections of Si, 2 along the total angular 
momentum axis, and the angle between S\ and S 2 . To 
completely specify the unit vectors along Si and S 2 , we 
require one more angle. One can arrive at a convenient 
choice for this angle by recalling that in the equal mass 
limit, the total spin simply precesses about the total an- 
gular momentum axis. So we choose the remaining angle 
to be the angle between the projection of S perpendicu- 
lar to J and the projection of So perpendicular to J. In 
the equal mass limit, this angle is simply given by ip [cf. 
Eq. (|4.6[) ]. We now determine an approximate formula for 
its equivalent in the general case. For unequal masses, 
the evolution equations for S and A can be shown to be 



dS_ 

# 

dA 

dip 



1 

2 L 

1 
2 



0+ J xS + 0-JxA + a-AxS 



(4.48a) 



0+J xA + (3-JxS + a+AxS 



(4.48b) 



where a± = a\ ± a 2 and 0± = 0\ ± 2 . We next take 
another derivative of the evolution equation for the total 
spin, and project the result perpendicular to J to obtain 



d 2 S A 
dip 2 
d 2 A x 



-uj^S 1 - = k s A a 



2 A _L c ± 

#2 , u A A = ^ A S , 
where uj 2 s A and k$,a are constants given by 



(4.50a) 
(4.50b) 



J S,A 

1 

4 



1X>" 



2 + +0l + (2a_0+ + a+0-){J ■ A) 



-a-(3-{j-S)+a 2 _{A 2 )-a-a+{Sl-Sl) , 

(4.51) 

0% + 2 _ - (2a+0+ + a-/3-)(J ■ S) 



+a+/3_(J • A) + a 2 + (S 2 ) - a_a+(,S? - Sf ) 



(4.52) 



20 + 0_ - {2a_0+ + a+0-)(J ■ S) 



(J • A) + a+a- (S 2 ) - o?_ (Sf - Sf ) 



(4.53) 



20+0. + (2a+0+ + «_/?_) (J • A) 



-a + /3_ (J ■ S) + a + a_ (A 2 ) - a\(S 2 - S 2 ) 



(4.54) 



The various averages appearing above are 



2 O-L 



d 2 S 



dip 2 



1 r 

'4 



0\ + 2 _ + (2/3 + a_ + a+0-)(J ■ A) 



-a-0-(J ■ S) - a_a + (S • A) + a 2 _A< 

-~ \^20+0- - (2/3+cl + a+0-)(J ■ S) 

+a-0-(J ■ A) - a 2 _(S- A) + a_a+S 2 

(4.49a) 

The equivalent equation for A^ is obtain by interchang- 
ing S <-> A and a± <-* — a T . The coefficients of S 1 - 
and A 1 - depend on time only through the function p(ip), 
since J • S^g and S*i ■ S2 are all directly proportional to 
p, which is an oscillatory function of ip. Our key approx- 
imation is therefore to replace the coefficients of S 1 - and 



(J-S) 
(J -A) 

(S 2 ) 
(A 2 ) 



{J'S) + a-(p), 
(J ■ A) + a + (p), 
Sl + 20_(p), 
W- (P), 



A 2 
^0 



(p) = SiS 2 



Ao-e- 



2w 2 



(4.55) 
(4.56) 
(4.57) 
(4.58) 

(4.59) 



7 Since the evolution of p is on the same timescale as the preces- 
sion frequencies of S 1 - and A- 1 -, one might object that this not a 
good approximation. However notice that the magnitude of the 
terms proportional to p in 1 14. 491 1 is small (of order ~ aS) com- 
pared to the dominant constant terms. Thus we believe that this 
should be a reasonable approximation, at least to obtain some 
trial solutions which can then be compared against numerical 
integrations. 
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By postulating harmonic solutions, one finds the normal 
frequencies u>± of system (I4.50P to be 



0J± 



1 

71 



1/2 



(4.60) 

The general solutions for and A x are thus lin- 
ear combinations of harmonic functions oscillating 
at frequencies lu±. In terms of initial conditions 
Sq, Aq,(1Sq /dtp, dA^/dip, the solution for S 1 - is 



s 1 - = 



(u 2 s - u)i)S^ - k s A^ cos($+) 



(lj 2 s - &l)S£- - k s A£ cos($_) 



,/c 

I 2 »2 \ ""-'0 ""0 



dip 

I 2 -2 \ 



dAj 



sin(<IL 



sin($_) 



(4.61) 



where the phases are 



Jo 



(4.62) 



Then the angle T between S ± and Sq is obtained di- 
rectly from (|4.6ip as follows 



tanT 



S 1 - • Sjf ' 



(4.63) 



We have now assembled all the pieces necessary to write 
down final solutions for each spins. These are 



Si. 



s 



1,2 



(cos 0i,2) J + (sin 0i,2 cos ^1,2)^0" 



+ (sin0i j2 sin <£i i2 ) J x S 



(4.64) 



The angles parametrizing the direction of unit vectors 
along each spin are the following 



0i = arccos 


(aiS 2 )u + (J ■ 5i)o 


1 


(4.65a) 


<Pi = - 


f <P-), 




(4.65b) 


02 = arccos 


- (a 2 S 1 )u + (J ■ S 2 ) 


, (4.65c) 


¥2 = ^(<y5+ - 


- <P-), 




(4.65d) 



where 

¥>+ 

<P- 
with 



= 2 arctan 



arctan 



Di tanT - D{ 
D[ + D 2 tanT 

—duj dip 



(4.66a) 



(Si ■ S 2 ) + (J3-)u - cos 0i cos 2 

(4.66b) 



£>i = (5isin0! +S , 2 sin0 2 ), 



Do 



(Si sin 0i — 52 sin 02) tan 



(4.67a) 
(4.67b) 



Let us now look at some features of this solution. First 
of all, since the evolution of the components of the spins 
is described by non-linear oscillators (either with cubic 
potentials or time-dependent frequencies), the true fre- 
quency spectrum in the case of unequal masses contains 
an infinite number of components. However our approx- 
imate solutions demonstrate that only a few of these 
components are sufficient to obtain very good agreement 
with solutions obtained numerically, as depicted below 
in Fig[TJ More precisely we included here five different 
frequencies into our solutions for the spins, namely the 
fundamental nutation frequency Q, p (d<&/dip), its first two 
harmonics 2fl p (d$>/dip) and 3fl p (d& / dip) , and the two 
frequencies appearing in S 1 - and A which are u)±Q p . 

To illustrate the accuracy of our approximate solution, 
we show below a plot of Si ■ S 2 and S 1 ' • Sq, the binary 
parameters being specified in the caption of FigQ] When 
specifying the initial binary configuration in the figure 
caption, the angles 0i.2 refer now to the angle between 
each initial spin and the initial orbital angular momen- 
tum, whereas the azimuthal angles ipi t2 describe the pro- 
jection of each initial spin in the initial orbital plane. 
They are not to be confused with the angles given in 
Eqs. (|4.65[) . which use the total angular momentum as 
the z-axis. By convention we always pick ipi = when 
specifying an initial binary configuration. We present 
Si ■ S 2 = cos 0i2 and S 1 - ■ Sq in Fig[T]to provide contrast 
with the equal mass case. In the equal mass limit, cos 0i2 
is a constant of motion (fluctuations over the orbital pe- 
riod are averaged out) and T increases proportional to 
time. Here clearly cos0i2 here evolves with time. Note 
also that the evolution of cos T contains very interest- 
ing structure, as it exhibits significant distortions from 
a pure sinusoidal shape. More specifically one can no- 
tice two turning points in the direction of precession of 
the component of the total spin perpendicular to total 
angular momentum. The first turning point occurs at 
~ 25 orbits and the second at ~ 30 orbits. Our ana- 
lytic solution captures this peculiar behavior rather well. 
We should point out however, that this is not a generic 
feature of the solutions. For most binary configurations 
we investigated, we did not notice turning points in the 
precession motion of the total spin. 
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-0.5- 



_ Cos 9j 


(numerical) 


. . Cos e 


(analytical) 


— Cos Y (numerical) 


■ ■ Cos Y (analytical) 







20 40 60 

Number of orbits 



80 



FIG. 1: We plot, as function of time in units of orbital pe- 
riod, the angles cos 812 = Si ■ §2 and cos T = S ± • So for a bi- 
nary with initial configuration M = IOMq, /orb = 20 Hz, xi — 
X2 = 0.8,6*1 = 6 2 = 30°, ipi = 0° and (p 2 = 90°. The mass 
ratio q = M1/M2 is here equal to 3/2. 



The agreement between our analytical prediction and 
numerical result for cos #12 is such that the numerical and 
analytical curves are indistinguishable on the figure. On 
this plot the maximum difference between the analyti- 
cal approximation and the numerical result for cos #12 is 
2.6 x 10~ 5 . This is evidence that our approximate so- 
lution to the cubic oscillator governing the evolution of 
cos #12 is very good. For cosT, one can distinguish the 
two curves by eye, but only barely. In this case, the 
maximum difference between the analytical approxima- 
tion and the numerical result is 6.32 x 10~ 2 . This is 
not quite as accurate as cos #12, but recall that the an- 
gle T is obtain by approximating the actual evolution 
equations, rather than developing an approximate solu- 
tion to the exact evolution equation. One could easily 
improve (|4.6ip by taking into account the leading oscil- 
lating terms in p in addition of considering its average 
value alone. This would generate terms oscillating at 
different frequencies than ui±, and complicate solution 
(|4.6ip considerably. For this first analytical investigation 
however, we believe the accuracy of (|4.6ip as it stands is 
sufficient. 

In Fig[TJ we track the evolution of the system for 100 
orbits at an orbital frequency of 20 Hz, so that it lies 
near the low frequency limit of the LIGO band. For a 
10 Mq binary, this is clearly unrealistic, as the binary is 
expected to evolve significantly due to radiation reaction 
on that timescale. Thus incorporating radiation reaction 
in our solutions is relevant, and we address this issue in 
the next section. 



V. EFFECTS OF RADIATION REACTION 

In the previous section we solved the precession equa- 
tions analytically, but neglecting radiation reaction. This 



assumption implied that the total angular momentum J 
and the quantity A are constants of the motion, allow- 
ing derivation of (approximate) analytic solutions. In 
realistic applications, the precession equations need to 
be solved over the radiation reaction timescale, and it 
is therefore important to discuss the corrections to the 
solutions previously obtained due to radiation reaction. 

In the presence of radiation reaction, the total angular 
momentum evolves according to 



dJ 

~dt 



32 pM 2 

-i(J-s] 



(5.1) 



In most astrophysical situations, the direction of J un- 
dergoes only a small precession around some average di- 
rection. This is known as simple precession Q. In rare 
cases it may happen that the total spin nearly cancels 
the orbital angular momentum, which may lead to wild 
changes in the direction of J. This is known as transi- 
tional precession Q . In this paper we will restrict atten- 
tion to simple precession alone and assume for simplicity 
that the direction of J remains constant throughout the 
inspiral. 

To illustrate our procedure to incorporate radiation- 
reaction into the solutions derived in the previous section, 
consider the precession of the total spin in the equal mass 
case. In absence of radiation reaction, recall the solution 
is given by 



S = [J ■ S ] J + S£ cos ip +{J x Sfr] sin V, (5.2) 



where 



^= / ^ 



1 

2d 3 



7-§* 



\J\df. 



(5.3) 



When J depends on time, the evolution equation for S is 
still given by (|4.5p . If the direction of J were exactly con- 
stant, then the solution would still be (|5.2p . However if 
J varies, then corrections terms must be added to (|5.2p . 
However in this paper, we shall implement radiation- 
reaction in the simplest possible way, which is given by 
the following prescription. One starts by writing down 
the spins for a given initial binary configuration using 
the solutions from section ITVl Then our solution which 
incorporates radiation-reaction is obtained by computing 
the phases ip,w and a in the equal mass case, or $ and 
<f>± in the general case by performing explicitly the inte- 
grals appearing in Eqs. (j43| . (j4T8|) and l[£T9]) . or (TQ0|) 
and (|4.62p respectively. All the coefficients of all trigono- 
metric functions appearing in the spins are assumed to 
be constants, determined by the initial conditions, i.e. 
we neglect possible evolution on the radiation reaction 
timescale of these coefficients. This approximation, along 
with the assumption that J remains constant, are the 
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main sources of errors in our solutions with radiation re- 
action included. However for a first analysis, we prefer 
to limit complications and leave improving the treatment 
of radiation-reaction beyond what is done here for future 
work. 

In order to calculate all required phases, the time evo- 
lution of the orbital frequency and of the magnitude of 
the total angular momentum are needed. We discuss 
these two topics next. 



A. Evolution of the orbital frequency 



We use the evolution equation for the orbital frequency 
accurate to 3.5 PN order, with spin effects included how- 
ever only up to 2PN order (see e.g. [8J), since this is 
the commonly used model employed in the data analysis 
community. It is given by 



96 , 743 - 9247? , , , n2/3 



12 ^ 



i=l,2 



\, (L \ -S;) j 113^ + 75/; 



4tt (Mw) 



/ 34103 13661 



V 18144 2016 



V + teV 2 ) {Moo) 



18 



U/3 



48 



X1X2 



[247( 



S 1 -S 2 )-721(L N -S 1 ){L N -S 2 ) (Mw) 4 / 



- — (4159 + 15876?i) (Mw) 5/3 
672 



/ 16447322263 1712 



16 



V 139708800 



541 
~896 ? 



5605 
2592' 



856 
105 



log 



16(Mw) 2/3 



7B H 7T 

105 ' 3 



/ 273811877 451 



V 1088640 



(Mw) 



\4032 



6048 



1512 



r 



91432\ 
13860 ) 



/4415 358675 91495 2 1 ,„ n7/3 , , . 

-r? - — — — — ry tt(Mw) 7/3 (5.4) 



where ^1,2 = |Si,2|/M 2 2 are the dimensionless spins, 
rj = fi/M is the symmetric mass ratio and 7^ is the 
Euler-Mascheroni constant. It is worth mentioning here 
that there are other contributions to (|5 .4[) from the spins 
at 2PN order, namely the quadrupole-monopole interac- 
tion given by Poisson [13] , and radiative self-spin contri- 
butions computed by Mikoczi, Vasiith and Gergely [23j |. 
There are also other contributions from spins at higher 
PN order which we also dropped here. Clearly it is for- 
mally inconsistent to neglect all these contributions, we 
nevertheless decided to do so for simplicity, and also to 
compare directly with a known model recently studied in 
the context of data analysis Q. This comparison should 
still be useful in estimating the errors generated by us- 
ing our analytic solutions in (|5.4p instead of the solution 
provided by numerical integration. 

Formulas for the total angular momentum and spin 
phases below involve the following function 



I = 



32 77 r 



(Mw) 8 / 3 



dui. 



(5.5) 
(5.6) 



Substituting (|5.4[) into (|5.6[) and performing a Taylor ex- 
pansion accurate to 3.5PN, one can in principle perform 
the integral, provided that some expression for the spins 
is provided. To simplify the problem, we use some sort 
of average values for the quantities hjq ■ Si i2 and S ■ Si t2 
to avoid introducing complicated non-linearities in the 



problem. More precisely, we replace in (|5.4[) Ln ■ S\$ 
by (J ■ §1,2) and Si ■ S 2 by (Si ■ S 2 ), each average value 
being computed from Eqs. (|4.55[) - (|4.59p . Our final result 
for the function / which appears throughout this section 



1 + H Cnl 



dy 

y 



Hy/y ) + J2^[y n -yS 



where 



Writing (|5.4p as 



y = (Mw) 1/3 . 



W 96 .r/o 

— = — t? Mw 5 / 3 
5 



y 



n=2 



(5.7) 



(5.8) 



(5.9) 



where the 6„'s are read off directly from right-hand side 
of (|5.4p . the coefficients c„ appearing in (|5.7p are 



c 2 
C3 
c 4 

C5 

C6 
C7 



b\ - 64, 



-62, 
-63, 
2 

2M 3 - 6b, 

+ 63 + 26264 - b 6 , 
-3b 2 2 b 3 + 2b 3 b 4 + 2b 2 b 5 



(5.10a) 
(5.10b) 
(5.10c) 
(5.10d) 
(5.10e) 
(5.10f) 
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Note that be, and thus eg actually depends logarithmi- 
cally on frequency. Here we neglect this frequency de- 
pendence and use b 6 = b§(t = 0) instead, since the errors 
generated by this simplification are smaller that the er- 
rors coming from averaging the terms depending on the 
spins in 63 and 64. 



B. Evolution of magnitude of total angular 
momentum 

We next move on to the computation of the magnitude 
\J\ of the total angular momentum. This is required in 
order to compute the phases ip, zu, a, $ and <I>±. 



1. Equal mass case 

In the equal mass case, the evolution equation for the 
total spin has the form 



dS 

—— oc J x S, 

at 

From (|5.1[) and (|5.1ip , we obtain directly 



(5.11) 



d{J ■ S) 

dt 



-7 [(J • S) S 2 } 



Equation (15.12p can be solved as 



J-S = er I [{J-S) -S 2 ] + S 2 , 



(5.12) 



(5.13) 



where / is given by Eq. (|5.7j) . Next contracting (|5.1| with 
J, one obtains 



2. General case 

We now repeat the computation for the case of arbi- 
trary masses and spins. Here however it is not possible 
to obtain J • S and hence J 2 exactly, since the evolu- 
tion equation for J • S contains a term proportional to 
J ■ (Si x S2), whose exact evolution is determined in 
terms of an elliptic integral. We start from 



djj ■ S) 

Jt 



- 7 [(J-S)-S 2 } 

+?(ai-a 2 )J-(SixS 2 ) 
at 



= -7 [(J ■ S)- S 2 ] + \J\( ai - a 2 
= - 1 [(J-S)-S 2 ] 



dp 
~dt 

3(M 2 - Mi) dp 



7// 



dt 

(5.16) 



where (|4.22p and (|4.26j) have been used, with terms of 
order A neglected in (|4.22|) . The solution to (|5.16p is 



'-^dt' 



J S = e~\j ■ S)o + e~ 7 I e'-fS 

Jo 



, 3(M 2 - Mi) _j f* jdp 

+ Tp e L e dv dt - (5 ' 17) 

Since p is a function oscillating over the precession 
timescale, we will assume that the remaining integral in 
(|5.17[) averages out to a negligible contribution, and that 
we may substitute S 2 by its initial average value (S 2 )o 
[cf. Eqs. (|4.57p and (|4.59|) ]. The remainder of the calcu- 
lation mirrors the manipulations of the previous subsec- 
tion, which lead to 



dJ 2 

dt 



= -2 7 [J 2 - J ■ S] , 



whose solution is directly given by 



-21 



-11 



Jo +2 



e 2I -J ■ Sduj 



1 Wo 

J 2 + 2(e T - 1) (J ■ S) c 



(5.14) 



ITS 2 



(5.15) 



where (|5 . 13[) has been used to obtain the second line. 
Note incidentally that Eq. (|5.15[) is exact, and thus does 
not rely on the assumption of constant J used throughout 
the paper. It is therefore also valid for spin configurations 
producing transitional precession. Since transitional pre- 
cession occurs when \ J\ — > 0, one can use (|5.15p to pre- 
dict if transitional precession can occur given an initial 
binary configuration (restricted for now to equal mass). 



J 2 = e- 21 



J 2 + 2 (e 1 - 1) (J • S) + {e 1 - I)' <S 2 ) 

(5.18) 



This expression can be used to predict (roughly) the pos- 
sible onset of transitional precession in the unequal mass 
case. 



C. Phases 

To complete the computation of radiation reaction ef- 
fects, we need expressions for the phases ip, tr, w, $ and 
$-t appearing in the solutions for the spins. 



1. Equal mass case 

The phases relevant to the equal mass solution are ip, w 
and a [cf. Eqs.flU]), (PTg|) and 1)435]) ] . We start with 
the phase ip, which is given by 
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2d 3 



2M 



\J\dt 



[4 + 2 (e> - 1) (J ■ S) + (e' - 1) 2 S 2 ] * 



doj 



(5.19) 



An exact evaluation of the above integral is impossible. 
It must therefore be evaluated numerically. We next turn 
our attention to the phase w, which is given by 



VJ 



o 2d 3 



\J\dt, 



n- 



2Af 



x [J 2 + 2 (e i - 1) (J ■ S) + {e 1 - l) 2 S 



2 c21 3 



(5.20) 



where 




0.025 0.05 0.075 

orbital frequency (in units of 1/M) 



0.1 



FIG. 2: This is a plot of the precession frequencies appearing 
in the spins, for the equal mass case, as function of orbital 
frequency. The binary configuration at orbital frequency of 
20 Hz is specified by the following parameters: M = 10 Mq , 
Xi = X2 = 0.8, 6i = 6 2 = 30°, fa = 0° and fa = 90°. The 
precession frequencies are expressed in units of instantaneous 
orbital frequency. The frequencies shown in the legend are 
defined as fii = dip/dt, fi2 = dzn/dt, £1$ — da/dt , fi4 = 
fii + CI3 and f7g = fii — ^3. 



n 



1 - 2a(J • S) + a 2 S 2 



1/2 



(5.21) 



Similarly to %j>, the phase w must be evaluated numeri- 
cally. Lastly we look at the phase <r, given by 



7-§> 



\J\dt 



Jo 2d 3 

r ssoj 2 u n duj 



(5.22) 



which we also evaluate numerically. To complete our 
analysis of the equal mass case, we give a plot of the 
precession frequencies as function of orbital frequency for 
the following initial binary configuration: M — 10 M Q , 
Xi = X2 = 0.5, 1 = 9 2 = 30°, fa = 90° and fa = 0°. 
Above M is the binary's total mass, xi,2 are the dimen- 
sionless spins and the angles #1,2 and fa t 2 are the polar 
and azimuthal angles giving the initial orientation of each 
spin in a coordinate system where the z axis is given by 
the initial orbital angular momentum. The zero of the az- 
imuthal coordinate (f> is defined by the initial component 
of Si lying in the initial orbital plane. 



2. General case 

The phases $ and $± appearing in the general case 
are defined in Eqs. (|4.40p and (14.62)) respectively. We first 
look at $, which is 



$ = 



w 
2d? 



5[w 2 xl + {x' ) 2 l 



\J\dt'. 
(5.23) 



Similarly to the phase w [cf. Eq. (|5.20[) ]. since w and e 
are now functions of frequency, one can transform the 
above integral into the following 



OJ 

2M 



w 



b[w 2 x 2 + (x' ) 2 ] 
6w 6 



7e~ 



i duj 



x [J 2 + 2 (e 1 - 1) (J • S) + (e 1 - 1) 2 S 2 ] 5 ™ 
L v ' J OJ 

(5.24) 

which needs to be evaluated numerically. The time de- 
pendence of w is obtained from (I4.44|) by substituting 
the time-dependent quantities 0^2 (J) and (3x^(1) into 
(|4.45p and (|4.47[) . The last elements required to con- 
struct our approximate analytic solutions for the spins 
are the phases $±. They are given by 



*+ = 



OJ 

2M 



- 7; A o 



oj± x 



1 doj 



[J 2 + 2 (e 1 - 1) (J • S) + (e 1 - l) 2 S 2 \ 5 — , 

0J 

(5.25) 

which cannot be evaluated in closed form. To complete 
our analysis of the general case, we give a plot of the 



g 0.4^ 




0.025 0.05 0.075 0.1 

orbital frequency (in units of 1/M) 



FIG. 3: This is a plot of the precession frequencies appearing 
in our approximate spins, a mass ratio q = M\jMi = 3/2, 
as function of orbital frequency. The binary configuration 
at orbital frequency of 20 Hz is specified by the following 
parameters: M = 10 M , xi = Xz = 0-8, 61=62 = 30°, $1 
u and (f>2 = 90°. The precession frequencies are expressed 
in units of instantaneous orbital frequency. The frequencies 
shown in the legend are defined as Qi = d&/dt, Q2 = 2d$/dt, 
SI3 = 3d<fr/dt , SI4 = d<fr+/dt and fis = d$-/dt. Note however 
that the true frequency spectrum contains an infinite number 
of components, but these five are sufficient to obtain good 
accuracy. 

precession frequencies as function of orbital frequency for 
the same binary configuration as the one for FigJ^l but 
now with a mass ratio q = M1/M2 = 3/2. 

This completes our construction of approximate solu- 
tions for the spins in presence of radiation reaction. We 
conclude this paper's analysis by comparing our analytic 
solutions to results obtained via numerical integration of 
the precession equations, including radiation reaction. 



D. Comparison with numerical solutions 

In this section we compare our approximate analytic 
solutions with solutions to the spin precession equations 
obtained numerically. We begin by simply plotting the 
components of the unit vector along Si, for the same 
binary configurations used to generate figures O and O 
to give a rough illustrative estimate of the accuracy of 
our approximate solutions. Figures 2] and [5] are gener- 
ated as follows. First, we integrate the orbital frequency 
evolution equation coupled with the precession equations 
numerically to generate the function u) num (t). Since our 
spins are viewed as functions of frequency when includ- 
ing radiation-reaction [cf. computation of the phases in 
previous subsection], the quantities that are plotted in 
Figs. H] and [S] are the components Si^Maum^)]- Thus 
we here first verify the accuracy of the functional form of 
our spins. 

Figures[3]and[5]both show excellent agreement between 
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FIG. 4: This is a plot of the components of the unit vector 
Si, for an equal mass binary. The binary configuration is 
the same as in Fig[2] The z-axis is the initial orbital angular 
momentum axis, and the x-axis is along the initial direction 
( 




10000 20000 30000 40000 50000 60000 
Time (in units of total mass) 



FIG. 5: This is a plot of the components of the unit vector 
along §1, for a mass ratio q = M1/M2 = 3/2. The binary 
configuration is the same as in Fig[3] Again the z-axis is the 
initial orbital angular momentum axis, and the :r-axis is along 
the initial direction of the component of Si lying in the initial 
orbital plane. 



the numerical results and analytical predictions. Over 
most of the inspiral, the relative difference between the 
numerical and analytical solutions is less than a percent, 
except near the end of integration, which we set when 
w rb = 0.1. Also one can notice in Fig. [5] that around 
t = 2.5 x 10 4 M, the discrepancy between the analyti- 
cal and numerical solutions gets a bit worse. This hap- 
pens when the total spin precession goes through turning 
points, as pointed out in the previous discussion of Fig[T] 
This effect is simply not captured as well as the rest of 
the evolution of the spins by our analytic solutions. One 
would need to incorporate more terms oscillating at dif- 
ferent frequencies when solving (|4.49[) in order to improve 
our solutions in that regime. Near the end of integration, 
the errors in the spin components are generically in the 
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range - 10% - 25%. 

We now move on to a different comparison, which is 
performed as follows. We first compute the orbital fre- 
quency as function of time by solving (|5.4p and the spin 
precession equations including radiation reaction simul- 
taneously, the result being a function of time denoted as 
w num (i), as before. One can then invert this equation and 
obtain time as function of frequency, i.e. t = t(uj nmn ). 

We now turn to coupling our analytic solutions for the 
spins with (15. 4)1 . To do this, we simply substitute the fre- 
quency dependent spins constructed in this section into 
the right-hand side of (|5.4|) , and then integrate the result- 
ing differential equation with respect to time numerically, 
the result being a function of time denoted as w an (<). We 
found generically that the coalescence time t an where w an 
diverges is shorter than the coalescence time associated 



with oj„ 



This will be relevant below. 



To compare these two solutions, we look at the dif- 
ference in accumulated orbital phase A<I> as function of 
frequency, defined as 
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q = 1:1 ; total cycles = 89.6 
q = 3:2 ; total cycles = 93.5 
q = 2:1 ; total cycles = 101.2 
q = 4:1 ; total cycles = 141.9 
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FIG. 6: This is a plot of accumulated orbital phase difference 
as function of orbital frequency for 4 different mass ratios 
(indicated in the legend) , and the same binary parameters as 
in FigH 



A$(w num ) = 



t("„um) 



(i') - uj an (t')] dt'. (5.26) 



We use the orbital phase for our comparison instead of 
gravitational wave phase here, since a typical waveform 
for a spinning black hole binary contains many different 
frequencies, which depend on the binary configuration. 
Thus to make our comparisons as uniform and simple as 
possible, we restrict attention in this paper to the orbital 
phase only. In Fig [6] we show a plot of A$/$ num as 

function of w n um, where $ n nm = J w n um(i') dt' ', for the 
binary configuration corresponding to Fig. [21 and for 
different mass ratios. The results for a few other points 
in parameter space are summarized in table HI the data 
of Fig[6] corresponding to the first block reported in the 
table. In Table [H the quantity A$/$ num is evaluated at 
the coalescence time t an . 

Generically, the errors are larger for bigger mass ratios 
since it takes more orbital cycles to sweep through a given 
frequency band. Thus there is more time for the errors 
inherent to our approximation scheme to build up. We 
believe this simple argument explains the generic features 
of Fig|6j and Table Q] The relative phase errors range 
from about 0.3% to about 1.8% for - 80 - 140 orbital 
cycles. This corresponds to errors of order ~ 6 radians 
over ~ 100 orbital cycles, which is a little too large for 
direct implementation in template banks, which require 
~ 1 radian of accuracy over ~ 500 orbital cycles. Thus 
our solutions would need additional refinement before be- 
ing suitable for template bank generation. However we 
believe that this ~ 1% accuracy is quite sufficient for a 
first cut at comparing the results of numerical relativity 
with post-Newtonian theory. 



TABLE I: This table displays the results of our comparison 
for a representative sample of binary parameter space. All 
binary configurations have a total mass M = 10 Mq and the 
initial spins have xi = X2 — X- The initial azimuthal angles 
are given by ipi — 0° and <^2 = 90°. 



binary 
parameters 



X = 0.8 

01 = 30° 

02 = 30° 



X = 0.8 

01 = 90° 

02 = 90° 



X = 0.8 

01 = 150° 

2 = 150° 



x = i.o 

01 = 30° 

02 = 30° 



x = i.o 

01 = 90° 

02 = 90° 



x = i.o 

01 = 150° 

02 = 150° 



mass ratio 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



1 : 1 

3 : 2 

2 : 1 

4 : 1 



number of 
orbital cycles 



89.6 
93.5 
101.2 
141.9 



84.7 
88.2 
95.3 
133.1 



79.7 
83.1 
89.8 
125.1 



90.9 
94.9 
102.8 
144.3 



84.7 
88.2 
95.4 
133.1 



78.5 
81.8 
88.5 
123.2 



10 2 



-0.42 
-0.30 
-0.41 
-0.67 



-0.85 
-0.85 
-1.38 
-1.58 



-0.39 
-0.46 
-0.54 
-0.89 



-0.53 
-0.38 
-0.52 
-0.86 



-1.15 
-1.08 
-1.35 
-1.84 



-0.55 
-0.62 
-0.76 
-0.95 
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VI. CONCLUSION 

We presented a detailed analysis ol the spin precession 
equations in black hole binaries, when the quadrupole- 
monopole contribution is taken into account. We showed 
that the (orbit-averaged) precession equations supple- 
mented by the quadrupole-monopole term possess a con- 
served quantity, which had not been previously noticed. 
The existence ol this conserved quantity allowed us to 
solve the precession equations exactly lor arbitrary spins 
in the equal mass case, neglecting radiation reaction. 
When masses are unequal, we resorted to a perturba- 
tive expansion to solve the precession equations approxi- 
mately, as an exact solution could not be obtained in that 
case. We then showed how to incorporate radiation reac- 
tion effects into our solutions for the spins adiabatically. 
To assess the accuracy of our approximate analytic solu- 
tions in the context of gravitational wave data analysis, 
we integrated the evolution equation for the orbital fre- 
quency accurate to 3.5PN order using our analytic spins, 
and also using the solution for the spins obtained from 



numerically integrating the precession equations. We 
showed that the relative errors in the accumulated or- 
bital phase is of order ~ 0.3% - 1.8% over - 80 - 140 
orbital cycles, the errors being the worst when the spins 
are initially lying in the orbital plane. Our solutions how- 
ever turn out to be quite complicated, which will prob- 
ably make them unattractive for direct implementation 
in template families. Nevertheless we believe that they 
are potentially quite useful for phcnomcnological analy- 
sis of numerical models of spinning binary black holes, 
e.g. those provided by the numerical relativity commu- 
nity, and to improve general understanding of precession 
phenomena in binary systems. 
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